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1 Summary 


A non-ideal gas model has been developed and retro- fitted into the MSES viscous/in viscid multi- 
element airfoil program. The specific applications targeted are compressible airfoil flows in wind 
tunnels employing heavy gases. The particular gas modeled in this work has been sulfur hexafluoride 
(5E6 )i although most heavy gases could be implemented if adequate state and caloric data were 
available. 

Numerical predictions with MSES indicate that the non-ideality of SF e significantly influences 
airfoil behavior in transonic flows, especially at the higher total pressures envisioned for pressurized 
tunnels. The dominant effect is that for a given freestream Mach number, local Mach numbers 
in supersonic zones are lower, and shocks are correspondingly weakened. Another (but apparently 
smcdler) effect is that for a given edge Mach number, a boundary layer in a heavy gas is theoretically 
somewhat more resistant to an adverse pressure gradient due to reduced adiabatic heating near the 
wall. 


As pointed out by Wagner and Schmidt [1], transonic small- disturbance theory is valid for 
non-ideal gases. Similarity between two flows can be obtained if the transonic similarity parameter 

K = Lll?— 

[M 2 ( 7'+l)] 2 / 3 

is matched, and if the pressure coefficients are scaled by the factor 

M 2 

A 3 


so that the quantity ACi must also be matched between the two flows. The parameters K and A 
above are defined in terms of an “equivalent” ratio of specific heats 7', which is derived in Appendix 
B for the second-order small- disturbance formulation employed in MSES. 

Although similarity between ideal and non-ideal inviscid transonic flows is rigorous in the con- 
text of transonic small-disturbance theory, a similarity rule cannot be formulated for viscous tran- 
sonic flows. In addition to the Reynolds number Re, Appendix C shows that an additional pa- 
rameter 7 V is introduced. This depends on the gas properties and local Mach number, and scales 
the effect of the local Mach number on the displacement thickness. It therefore affects viscous 
displacement effects and boundary layer response to pressure gradients in compressible flows. It is 
highly unlikely that the parameters M , 7', Re, and 7 V can all be combined into one similarity rule 
for viscous transonic flows. Fortunately, numerical experiments indicate that matching if, ACi , 
and Re (or Af*, Cl, and Re) still gives good correspondence between air and heavy-gas flows. 
Apparently, the effect of 7 V is not nearly as significant as the other three parameters. 


Figure 1 compares C p vs xjc curves for the RAE 2822 airfoil [2] at M - 0.735 for air, for 5F 6 
at 1 atm, and for SFe at 3 atm. Figure 2 makes the comparison at a fixed M * = 0.765 instead 
of a fixed M . Figure 3 in turn makes the comparison at a fixed K and ACl (corresponding to 
M — 0.735 and Cl — 0.743 for air). Clearly, matching M* or K is more appropriate for evaluating 
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Figure 1: C p distributions for RAE 2822 airfoil at M - 0.735 for air, SF e (1 atm), and SF 6 (3 
atm). Ci — 0.743, Re — 6.2 million. 

transonic flow characteristics. To illustrate further, drag- divergence behavior for air, SFe (1 atm), 
and SFe (3 atm) is shown versus Af and M* in Figures 4 and 5. As expected from the C p 
comparisons, the effects of the type of gas on transonic drag rise are much smaller if Af* is used as 
the compressibility parameter in lieu of Af. The Mach sweep results were not performed at fixed 
K and ACi , since it is not clear how to scale the profile drag coefficient Co over this sweep. In 
principle, the pressure drag should be seeded by A, while the friction drag should perhaps be left 
unsealed. However, it is impossible to separate these drag components in an experiment, since only 
the total drag is obtained from a wake survey. 

For high-lift configurations, small-disturbance theory is obviously invalid, but numerical studies 
indicate that matching K and ACi (or alternatively matching M* and Cl) still gives a reasonably 
good match between air and heavy-gas flows. Figure 6 shows the inviscid C p distributions over a 
slatted two-element airfoil described in reference [3]. A freestream Mach number of Af = 0.30 in 
air produces a fairly strong shock on the slat find a somewhat weaker shock on the main element. 
Figure 7 compares the C p distributions on the slat for the three gas cases at a fixed sonic Mach 
number Af* = 0.3257 (corresponding to Af = 0.30 for air) and Cl = 2.85. The comparison is quite 
reasonable. It should be stressed again that simply matching the usual freestream Mach number 
Af = /a x and unsealed Cl gives a very poor match in all cases, except of course in effectively 
incompressible flows where any gas non-ideality is irrelevant. 


2 


3 














The bulk of the heavy-gas model development and application to transonic, inviscid flows is 
documented in the SM Thesis of Marc Schafer, which is attached as Appendix A. As mentioned 
previously, Appendix B derives the farfield behavior of a non- ideal airfoil flow. This was required 
for implementation of new outer boundary conditions for the MSES code. Appendix C derives 
the shape parameter compressibility correction for an adiabatic boundary layer in non-ideal flow. 

This was required to implement new heavy-gas correlations for the MSES integral boundary layer 
formulation. 
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Appendix A 



Modeling of Heavy Gas Effects on Airfoil Flows 


by 

Marc Alan Schafer 


Submitted to the Department of Aeronautics and Astronautics 

on May 3, 1992 

in partial fulfillment of the requirements for the degree of 
Master of Science in Aeronautics and Astronautics 


Thermodynamic models were constructed for a calorically imperfect gas and for 
a non-ideal gas. These were incorporated into a quasi one dimensional flow solver to 
develop an understanding of the differences in flow behavior between the new models and 
the perfect gas model. The models were also incorporated into a two dimensional flow 
solver to mvestigate their effects on transonic airfoil flows. Specifically, the calculations 
simulated airfoil testing in a proposed high Reynolds number heavy-gas test facility. The 
results indicated that the non-idealities caused significant differences in the flow field, 

but that matching of an appropriate non-dimensional parameter led to flows similar to 
those in air. 
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Associate Professor of Aeronautics and Astronautics 



4 


0 



Acknowledgments 


would bit. to express my thanks all thoso .ho mad. .his thesis possibU. Fist, 
to Mask Dreia whose brilliance and ingenuity have _,d as an inspiration in ^ of 
my studies. Also, to Harold 'Guppy' Youngran -hose laadarship during the Daadalu, 
project helped me to realise what it really means to be an engineer. 

I would also like to thank my parents and the res, „f m , fa^y. WitW ^ 
support, I never would have made it as far as I have. 

My appreciation also g„„ to the NASA Ames research center and the NDSEG 

e owshtp program without whose financial support this thesis would never have hap- 
pened. y 


PRECEDING PAGE BLANK NOT FILMED 


5 




6 



Contents 


Abstract 

Acknowledgments 

1 Introduction 

2 Real Gases 

2.1 Calorically Imperfect Gases . . 

2.2 Non- Ideal Gases 

3 Solving the Euler Equations 

3.1 Calorically Imperfect Gas . . . . 

3.2 Non-Ideal Gas 

4 Results 

4.1 One Dimensional Duct Flow . . . 

4.2 Two Dimensional Results .... 

5 Conclusions 

A Curve Fit For SF e State Equation 


3 

5 

10 

11 

11 

12 

15 

16 
17 

20 

21 

23 

26 

27 


PRECEDING PAGE ENJtf-JX NOT FILMED 7 


1 / 



B MSES Subroutine for Non-Ideal Gas Model 


28 


Bibliography 


43 



8 



List of Figures 


4.1 Stagnation Pressure Ratio( Strength) vs. Upstream Mach No. for Air 



and SF$ at latm and 3atm ..... 


90 

4.2 

One dimensional Duct Flow . . . 


. L\j 

91 

4.3 

Shock Strength and Location vs. /3 . . . 


• Z 1 

99 

4.4 

Shock Strength and Location vs. Zn . . . 


99 

4.5 

Comparison of Air and SF& at Fixed M and Ci . 

. 24 

4.6 

Comparison of SF 6 at latm and 3atm to Air, 

M * = .740, C L = . 9 . . 

. 24 

4.7 

Comparison of SF& at latm and 3atm to Air, 

AT* = .732, C L = .75 . 

. 25 

4.8 

Comparison of SF 6 at latm and 3atm to Air, 

k = .439, AC l = 2.18 . 

. 25 


9 



Chapter 1 

Introduction 


In the past few decades, the design and development of large transport aircraft has 
relied on wind tunnel data taken at significantly lower Reynolds numbers than those 
found in operation. The drawbacks of this subscale data become apparent when one 
considers phenomena such as attachment line transition or similar aspects of boundary 
layer behavior at high Reynolds numbers. 

The need for accurate wind tunnel data clearly mandates the construction of a 
suitable high Reynolds number test facility. However, the cost of building a large at- 
mospheric tunnel and large tunnel models is prohibitive. Higher Reynolds numbers are 
often achieved by pressurizing tunnels to effectively increase the density of the air. This 
alternative is practiced only up to a point. 

A potential solution following the same basic idea relies upon the use of gases with 
significantly higher molecular weights than air. Candidate gases include Freon- 12 or 
Sulfur Hexaflouride (SF 6 ), but the use of non-breathable gases clearly causes some 
problems. These problems will likely be insignificant to the cost and operational ad- 
vantages of such a facility. Combining heavy gases with pressurization would allow test 
Reynolds numbers comparable to those on large transports in flight [1]. 

One complication is that Freon and 5 F* have significantly different thermodynamic 
properties than air, especially at elevated pressures. Heavy gases do not follow the ideal 
equation of state P = pRT nearly as well as air does, nor do they maintain a constant 
ratio of specific heats 7 = Cp/c, over any significant temperature range. The following 
discussion will attempt to quantify the potential importance of these effects through a 
computational study. 



Chapter 2 

Real Gases 


The thermodynamic relations specifically subject to real gas effects are the state equa- 
tion 

P = P RT (2.1) 

and the caloric equation, 

h = J Cp dT = c p T (2.2) 

these particular forms only being valid for a perfect gas. Real gas effects may be divided 
into two categories: 

1. Calorically imperfect gases for which c p depends on temperature, but which still 
satisfy equation (2.1). 

2. Non-ideal gases for which c p depends on both pressure and temperature, and 
equation (2.1) no longer holds. 

The first effect results from the introduction of multiple vibrational modes for poly- 
atomic molecules which become more important at higher temperatures. The second 
effect depends on intermolecular forces which become stronger as a gas moves towards 
liquefaction, ie. higher pressures and lower temperatures. 


2.1 Calorically Imperfect Gases 


The only difference between a perfect and an imperfect gas stems from the dependence 
of c,, on temperature in the imperfect case. A cursory examination of experimental data 
for 5Fj shows that, in the range of temperatures likely to be found in a wind t unn el 



test, this dependence is linear in temperature. 


c p (r) = a + bT 

Therefore, equation (2.2) becomes 

bT 2 

h(T) = aT + 

which may be easily inverted to find T(h). 

- -5 + /(i) ! + T 

2.2 Non-Ideal Gases 

The state equation for a perfect gas (2.1) derives from a kinetic model of gas molecules 
which assumes that the molecules are point masses and that they do not exert any forces 
on one another except instantaneously during collisions. Clearly these assumptions 

become less accurate as the molecular weight of the gas increases. Van der Waals’s 
equation 

(p + p 2 aj (1 - p/3) = pRT (2.6) 

contains two correction to equation (2.1): a corrects the pressure to account for inter- 
molecular attraction, and (3 corrects for the volume of the molecules themselves. 

Using a non-ideal state equation like Van der Waals’s causes many serious compli- 
cations as enthalpy, Cp, 7, etc. now depend on pressure as well as temperature. Despite 
these complications, enthalpy and entropy must remain state variables regardless of the 
form of the state equation. That is, local entropy and enthalpy must depend only on 

the local pressure and temperature and not on the upstream conditions (ie. the gas 
history). 

Liepmann and Roshko [2] equate this condition with the requirement that a canonical 
equation of state must have one of these four for ms : 

e = «(*,P) (2.7) 


(2.3) 


(2.4) 


(2.5) 
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h = A(j,p) 

f = f(T,p) 
9 = g(r,p ) 


( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 


Here e _ h - p/p is the usual internal energy, / = e - Ta is the free energy, and 
g = h — Ts is the free enthalpy. 


For a conventional flow solver, the enthalpy defintion (2.8) appears best; however, 
specifying the state in this specific form is not convenient because the entropy s is not 
readily available to the flow solver. Liepmann and Roshko propose a more suitable form 

pRT = Z ( p,T ) (2.11) 

which requires T(p, h) to have a form which makes h a state variable. 


For a Van der Waals’s gas 


Z = — 1 -fg. 

1-/3 p RT 

which clearly approaches the ideal state equation for a ,/ 3 
values of a and (3 


( 2 . 12 ) 

0. For typically small 


(»3) 

where the second approximation is made to make Z = Z(p , r) explicitly. Liepmann and 
Roshko write equation (2.13) in more general form as 

z = < 2 -“> 

with p e and T c being the critical pressure and temperature of the gas, and <j> evidently 
being a universal function which they tabulate for gases other than air but with ap- 
proximately the same molecular weight. For heavier gases such as SF t it is best to fit a 

curve to experimental data as explained in Appendix A. For SF«, a good curve fit takes 
the form 

^(^f) = CJ +Ci(y) + c ° ( 2 - 15 ) 

It is now necessary to determine the specific heat capacity Cp(p, r) so that the enthalpy 
function h(p, t ) can be obtained. Liepmann and Roshko combine two forms of the 
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equation of state h(p, r) and a(p,r) into the fundamental reciprocity relation between 
h(p, t) and p(p , r) 


M = 1 _ r W£) 

dp P dT 


(2.16) 


which is valid for any gas. Combining this with the state equation (2.11) gives 


dh 

dp 


RT 2 


( dZ\ RT Ci .(t c \ 

\af), = T* (t) s < 217 > 

Since dhjdp = T{t) only depends on the temperature, both h and c p must be linear in 
the pressure as follows. 


Mp> t ) = f Cp(r) dT + pT(r) 

, \ dh 

C P (P,T) = — 

v , dT 

= C " (r) + P df 


(2.18) 

(2.19) 

( 2 . 20 ) 
( 2 . 21 ) 


As in the case of the calorically imperfect gas, Cp(r) has the form 


c "p(r) = a + bT (2.22) 

Substituting this into the enthalpy equation gives 

Hp, t) = aT + SI + #'(£i) (2.23) 

It is also possible to determine the caloric equation by expressing the internal energy 
(e) as e(p, t) [3]. 
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Chapter 3 

Solving the Euler Equations 


These gas models may be readily integrated into an existing Sow solver which solves 
the integral form of the steady Euler equations: 


/ 


pu- hdA 



u- hu + ph) dA 


ho S h + ^i 
2 


0 

0 

constant 


( 3 . 1 ) 

( 3 . 2 ) 

( 3 . 3 ) 


These equations are exact for any fluid flow, but must be supplied with a state equation 
to relate the pressure p to the enthalpy h and the density p. In addition, the upwinding 
scheme used to capture the shocks requires the local Mach number while the boundary 
conditions and evaluation of shock losses require the local stagnation conditions. 


It is desirable to nondimensionalize the equations, and the following scheme is used 
where () denotes the dimensional quantitiy and () re/ denotes a reference quantity: 

P ~ P/ftrf 
P = P/pTfd 
T = T/T^t 
h = 

Prttf 

Furthermore, Cp, c„, and R are nondimensionalized using R resulting in several new 
nondimensional par ame ters 


at = a/R 

0 = ^ 
2a 

» = Prtf/Pc 
r = Trrf/Te 
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For the results presented here, the reference conditions are chosen to be stagnation 
conditions. 


3.1 Calorically Imperfect Gas 


The nondimensional form of the caloric equation which governs the behavior of the 
imperfect gas is: 


h(T) = j Cp dT 

= aT + a(3T 2 


(3.4) 

(3.5) 


which may be inverted 


to give T as a function of h. 

T(h) = + y/ 1 + 4 0h/ 

K 1 2p 


a 


(3.6) 


With T obtained from h, p may be determined using the ideal gas law (2.1) and a 
specified value of p. The local Mach number comes from the familiar deftntion of the 
speed of sound: 


a 2 = dp 


dp 


= jT 


The local value of 7 may be found from equation (2.3). 


(3.7) 


Cp _ a + 2a(3T 
c. 1 o — 2a/3T 


(3.8) 


The last re ma i nin g difficulty is the determination of the isentropic relations between 
pressure, density, and temperature. These relations are necessary to calculate stagnation 
conditions from flow conditions. The familiar perfect gas relations 


T_ 

T 0 





P_ 

Po 


do not hold for a calorically imperfect gas. 




P_ 

Po 


(£)* 


The proper forms are obtained from the formal statement, 


dh = Tds 

P 


(3.9) 


HD 
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and for an isentropic process ds — 0: 


dh = ± 
P 


From the definition of enthalpy dh = c p dT, and for an ideal gas p/p = T, 
(3.10) becomes 

c p (r) dT dp 
T ~ 7 

Integrating this equation gives 


(3.10) 


so equation 


(3-11) 


7 = ei P(- Q logI’ + 2a^(l-r)) (3.12) 

and the isentropic density relation then follows directly from the state equation. 


P_ = P Tjho) 
Po Po T(h) 


(3.13) 


Strictly speaking, solution „f the Euler eqlJation! , tquir „ Dothillg ^ However 

a Newton- Raphson technique is nsed, all of the necessary equations must be linearized 

for the Jacobian matrix. In the case of the odorically imperfect gas, the equations 

are slightly more complicated than for a perfect gas, but they may still all be written 

explicitly. Therefore the linearizations are easily done by differentiating the relevant 
equations. 


3.2 Non-Ideal Gas 


The nondimensional eqtmtion. describing the non-id«d gas are the state equation 

JL _ l + P*4>(Jr>) 

pT Z 0 (3-14) 

and the caloric equation. 


A (P. T ) = \aT + a/3T 2 + )1 -L 

L r tt'\ Zo 

Z 0 is another parameter which may be described in terms of x and r. 

Z0 = ^ = 1 + »^) 


(3.15) 


(3.16) 
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The non-ideal gas presents some difficulty as the enthalpy depends on the temper- 
ature and the pressure. Therefore, from equations (3.14) and (3.15), p and T may be 
found using a Newton- Raphson system to drive the following residuals to zero. 


Ri(p, T) 

R 2 (p,t) 


p_ _ l + px4>(tt) 

pT Z 0 


h - 


aT + a/3T 2 + p-tff 

T 



Zo 


(3.17) 

(3.18) 


The local Mach number depends on the speed of sound which must be found from 
the definition: 


This is calculated as follows: 



(3.19) 


dp = d -l 
d P 


, dp 
dp + — 

, dh 


dh 


but dh = dp/ p for an isentropic process, and hence 


a 2 = dp 


dp 


1 ~ ^ I,* 


The local 7 really has no meaning and need not be calculated. 


(3.20) 


(3.21) 


The extra complexity of the non-ideal gas appears in the calculation of the sensitiv- 
ities. Since p and T are found by an iterative process they must be found by perturbing 
the Jacobian matrix of the converged Newton-Raphson system. A perturbation in h 
and p is related to a perturbation in p and T by the condition that the R(p, r, h, p) must 
remain zero. 


U.l 

► = 0 = 

9Rt 9Rt 

bh &p 

i 

[m) + 

^ 1& 

| Sp 

1 6R7 J 


§Rj SR, 

dh dp 


l ‘o] 


i «) 


Numerically inverting this system gives the required derivatives. 


U] 

j 


‘in 


i 

\ Sh 

ST J 





l * . 


The second derivatives are found in a similar fashion starting instead with 
as the residuals. Using a subscript notation for the derivatives (§g = p h ): 


(3.22) 


(3.23) 




(3.24) 
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A similar system with R Xp and R 2p as residuals is also formed. As above, numerically 

inverting gives fj£ = etc. These manipulations are implemented in 

the source code in Appendix B. 


The last remaining task is calculation of the stagnation conditions and, again, it 
is not possible to find an analytic expression. Another Newton-Raphson system is 
constructed where the first residual comes from equation(3.15): 


R\ = h Q - h(p,T) 

The second residual is derived by rearranging equation(3.9) 

, dh dp 

da = — -f — £- 

T pT 

= %dT + 

T T p 


= $iT + rf(pr-L^) -rdf,*)- ± 


Integrating gives: 


■’(P) T ) - J -jr dT + pK ^ j - ln(p) 

The second residual may then be formed 


( 3 . 25 ) 

( 3 . 26 ) 

( 3 . 27 ) 

( 3 . 28 ) 

( 3 . 29 ) 


.R 2 = Si - s(p, r) (3.30) 

where is the entropy of the static conditions. 


Driving these two residuals to zero gives the stagnation conditions po, T 0 . The 
derivatives fg*, etc, needed for the Newton-Raphson solver may then be found by 
perturbing the converged Jacobian matrix and relating the resulting derivatives to the 
static conditions through the chain rule and equations (3.15) and (3.29). This process 
is identical to the one used above to find p and T and their derivatives. 
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Chapter 4 

Results 


After developing the models for the calorically imperfect and non-ideal gases, the next 
step was to evaluate the differences these changes caused in inviscid flows. The primary 
quantities of interest are the location of shocks and their strength which is defined as the 
ratio of of stagnation pressures across the shock. For a perfect gas, the shock strength 
may be expressed as a function of the upstream Mach number M\. 


P02 _ 

1 + - 1) 

-l/(7-l) 

h + i)M? 1 

P01 ~ 

7 + 1 J 


[(7 -1)M* + 2 


However, for the non-ideal gas, this relation must be calculated numerically. 



Figure 4.1: Stagnation Pressure Ratio( Strength) vs. Upstream Mach No. for Air and 
SFs at latm and 3atm 
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Figure 4.2: One dimensional Duct Flow 


4.1 One Dimensional Duct Flow 


The first comparison of the different gas models was a study of the flow in a converg- 
ing/diverging nozzle using a quasi one dimensional Euler solver. This flow is character- 
ized by some flow at the throat with a shock downstream to match the specified exit 
pressure as shown in figuxe( 4 . 2 ). 

As a basis for comparison of the different gas models in a duct flow, the non- 
dimensional reference enthalpy (Aopo/po) was made equal for all three cases. 

*° = ( 4 -2) 
= a(l + /3) ( 4 . 3 ) 

a(l+ 0 ) + *^(l) 

Zo ^ ( 4 - 4 ) 

With ho held constant, 7 therefore depends on a, / 3 , t, and r. The exit presure ratio is 
also held constant. Under these conditions, the slope of the c, versus T curve (/ 3 ) had 

little or no effect on shock strength or position relative to the perfect gas as shown in 
figure(4.3). 
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Figure 4.3: Shock Strength and Location vs. (3 


For the non-ideal gas, it and r are not really independent parameters and may be 
combined into Z 0 . Figure(4.4) shows the variation in shock strength and position as 
functions of Z 0 and the corresponding perfect gas results with 7 adjusted to preserve the 
stagnation enthalpy as above. These plots clearly show that it is not possible to mimic 
the effects of the non-ideality by changing 7 as in the case of the caloricaUy imperfect 

gas. The difference in shock strength and position becomes larger and larger as the gas 
becomes less ideal. 



Figure 4.4: Shock Strength and Location vs. Z 0 
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The last test conducted with the one dimensional flow model was to determine the 
effects of the various gas models on the upwinding scheme needed for stability of the 
numerical scheme. The flow solver drives the momentum equation residual to zero, 

fli = PWiAiiqi - 9.-1) + Pi Ai - Pi-Ui.x + — + o A ~ 1 (Ai - A^) (4.5) 

where the upwinded speed is defined as 


* = ft - M<h - q<-i ) 


and ^ is non-zero only if M, is greater than M c . 


MM(ft)) = ^ 
7 



(4.6) 


(4.7) 


Initially, the exact 7 was calculated at each node along with all the necessary lin- ’ 
earizations and used in the upwinding scheme. Under these conditions, the flow solver 
converged with M c < 1 . However, the upwinding is relatively insensitive to the exact 
value of 7 even though the stability analysis used to derive equation(4.7) ignored 7 
perturbations. Using a constant value of 7 had absolutely no effect on the viable range 
for M e or the rate of convergence. 


4.2 Two Dimensional Results 


The subroutine which appears in Appendix B was incorporated into MSES, the multi- 
element version of the two dimensional transonic airfoil design/analysis code ISES [4]. 
Numerical experiments carried out were limited to single-element in viscid cases to more 
clearly demonstrate the effect of the new gas model. Figure(4.5) shows an overlay of the 
Mach distributions for a test airfoil run in SF t at two different stagnation conditions and 
in air. All three cases are at matched freestream Mach number and lift coefficient. Note 
that they are not at the same angle of attack. The SF 6 is characterised by stagnation 
pressures of latm and 3atm and a stagnation temperature of 310AT. 

Airfoils tests in heavy gases will be much more worthwile if some relationship may 
be found so that the tests reflect the airfoil performance in air. The only parameters 




Figure 4.5: Comparison of Air and SF e at Filed M and C L 

which may be adjusted in a wind tunnel test are the Mach number, stagnation con- 
ditions, and angle of atttack or C L . Figure(4.5) shows an attempted match keeping 
M and Cl constant: clearly, this is not an effective technique. After a good deal of 
experimentatation, the best match was achieved by running the different gases at the 
same AT which is defined as the ratio of freestream velocity to the speed of sound at 
some conditions. Figure(4.6) shows the case in air from figure(4.5) compared with SF 6 
(latm and 3atm) at the same M*. 



Figure 4.6: Comparison of 5 Fa at latm and 3atm to Air, M* = .740, Cl = .9 
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Figure 4.7: Comparison of SFq at latm and 3atm to Air, M* = .732, Cl — .75 

A case with a weaker shock, figure(4.7) was used to further verify this relationship. 
The match is slightly worse, but this is to be expected because a weak shock is much 
more sensitive to small changes in M than a strong one. As an alternative to matching 

M\ Anderson [5] proposes matching the small disturbance similarity parameter k and 
ACl where 

1-A& 

(Af£( 7 ' + 1 )) 2/3 (4.8) 

Mlb' + 1 ) 

1 - ( 4 - 9 ) 




Figure 4.8: Comparison of SF& at latm and 3atm to Air, k = .439, ACl = 2.18 
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Chapter 5 

Conclusions 


The models derived above adequately describe the thermodynamic behavior of non-ideal 
and calorically imperfect gases. Despite some minor complications in linearizing these 
models, they were implemented in routines suitable for incorporation into existing flow 
solvers based on Newton’s method. First, a quasi one- dimensional flow solver was used 
to examine the influence of the various non-dimensional parameters which govern the 
behavior of the different gases. 

Transonic airfoil test cases for air and SF 6 were then used to study the influence of 
parameters which may be controlled in a wind tunnel experiment: stagnation pressure, 
freestream Mach number, and angle of attack. The goal of this study was determine 
the conditions under which a wind tunnel test in a heavy gas would produce results 
comparable to those found in air. Matching Af* and C L or k and AC L were both 
effective for the test cases presented here. Further study is necessary to determine 
which is best for multi-element cases. 

The results are encouraging in that they definitely hint at the possibility of directly 
relating heavy gas test data to performance in air. It is first necessary to verify experi- 
mentally the model for SF 6 , and to investigate the effects of non-ideal gases on viscous 
flows. 



Appendix A 

Curve Fit For SF 6 State Equation 


A curve fit may be found for the function <f> (^) for any gas given experimental state 
data. With the density ( p ) measured at a number of different pressures (p) and tem- 
peratures (T), a vector is defined containing the difference between the real gas and a 
perfect gas at each data point. 


Z = 


At -1 


L eJilTr* - 1 


Defining 9 = , the matrix A contains the state information. 

Pi*? Pi#?' 1 ... Pi#\ -PxO i ^ 


A = 


PT'f 


(A.l) 


(A-2) 


Pm#m Pm # Hi 1 • ■ ■ Pm#m Pm#m Pm 

The goal is to find a state equation agreeing closely with the experimental data in 2 
but of the simple form: 


z(p,t) = 1 + -L- 

JVc/ L 


C n C n _i 


0 " 

ffn-l 


(A.3) 


Therefore 

2 ~ AC 

and (? is found by the technique of linear regression: 

<? = (A T A)~ 1 A T Z 


(A.4) 


(A.5) 


The results presented in this thesis were based on a quadratic fit for <(> from approximate 
data for SF&. The required data may be found in [6]. 
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Appendix B 

MSES Subroutine for Non-Ideal Gas Model 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 


c 

c- 


c 


c 


c 

c 

c- 


c 


subroutine hgparm(alfl ,btal , taul, cc0,ccl,cc2, hO) 


Initializes non-ideal gas routines. 

Formulation derived in Schafer SM thesis. 

Input : 

alfl Constants for Cp(T) in caloric equation: Cp = a(l + bT) 

betl 

taul Constant in phi(T) in non-ideality factor Z(p,T) 

ccO Constants defining phi(T) in polynomial form: 
ccl 

cc2 phi = cO + cl (tau/T) + c2(tau/T)**2 

Output : 

hO Enthalpy at reference conditions pO, TO 
Internal output: 

20 Ion-ideality factor Z(pO,TO) at reference conditions 


implicit real*4 (a-h,m,o-z) 
common /nongas/ 

* alf, bta, pi, tau, zO 
common /nonfit/ 

* c2, cl, cO 

put input parameters into common blocks 
alf = alfl 
bta = btal 

tau * taul 

cO * ccO 
cl = ccl 
c2 = cc2 

pi = 1.0 

calculate reference non-ideality factor and enthalpy 

zO s i.o + pi*(c2/tau**2 + cl/tau + cO) 

hO = (alf*(l . + bta) + pi/tau*phid( 1 . /tan) ) / zO 


: 3 ^ 
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return 

end 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 


subroutine nideal(hO ,r ,q, p , p _ r ,p_q, 

* msq,msq_r ,msq_q) 

Calculates pressure and Mach number for specified 
stagnation enthalpy, density, and speed. 


Input : 

hO stagnation enthalpy 
r density 

q speed 


Output : 
P 

P- r 

p-q 

msq 

msq.r 

msq_q 


pressure 

dp/dr 

dp/dq 

square of Mach number M"2 

dM~2/dr 

dM~2/dq 


implicit real*4 (a-h,m,o-z) 


c set static enthalpy 
h = hO - O.S*q**2 
h_q = -q 

c 


c sat prassuxa and tamparatura and darivativas 
call ngaspt (h , r , p , p_r , p_h , p_rr , p_hh , p_rh , 

* t ,t_r,t_h,t_rr,t_hh,t_rh) 

P-q = P_h*h_q 


c sat spaad of sound squaxad: a*2 = dp/dr (at constant s) 

•*q = p.r / (l. - p_h/r) 

asq_r = p_rr / (1. - p_h/r) 

* “ P- r / (1. - p_h/r)**2 *(p_h/r**2 - p rh/r) 

*«q_h = p_rh / (1. - p_h/r) 

* + P_r / (1. - p_h/r)**2 *p_hh/r 
asq_q = asq_h*h_q 

c 

set Mach number squared 
msq = q**2/asq 
msq.r = -msq/asq * asq.r 

®*q-q = -»sq/asq • asq.q + 2.*q/asq 

c 

return 

end 


subroutine 


ngaspt (h,r,p,p_r ,p_h,p_rr,p_tah,p_rh, 
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t ,t_h, t_rr ,t_hh,t_rh) 

Calculates pressure and temperature for 
specified static enthalpy and density. 


c 

Input : 


c 

h 

enthalpy 

c 

r 

density 

c 



c 

Output : 


c 

P 

pressure 

c 

P- r 

dp/dr 

c 

P-k 

dp/dh 

c 

P- rr 

d~2p/dr~2 

c 

p_hh 

d~2p/dh~2 

c 

p_rh 

d A 2p/drdh 

c 

t 

temperature 

c 

t_r 

dt/dr . . . etc 


c 


implicit real*4 (a-h,m,o-z) 

dimenaion a(2,2), ai(2,2), aih(2,2), air(2,2). 

* b(2,2), bh(2,2) , br(2,2) 
common /nongas/ 

* alf, bta, pi, tau, zO 

c 

c Newton convergence tolerance 

data eps /5.0E-6/ 

c 

c initial guess from imperfect ideal gas 
if (bta. eq, 0.0) then 
t = h/alf 

else 

t = (-1.0 + «qrt (1.0 + 4.0*bta*h/alf )) / (2.0*bta) 
andif 
p = r*t 


c lawton loop to convarga on corTact p,t 
itcon = 16 

do 100 itar=l, itcon 

c 

c aat and linaariza non-idaality factor Z(p,t) 
ttc = l./(tau*t) 
ttc_t = -l./(tau*t**2) 

c 

z = 1. + p*pi*phi(ttc) 

Z -P = pi*phi(ttc) 

z - t = p*pi*phid(ttc)*ttc_t 

c 

c raaidual 1: atata aquation 

raal = p/(r*t) - * / Z 0 

rl -P = l./(r*t) - z_p/z0 

rl_t = -p/(r*t**2) - z_t/z0 
c 
c 


'3Y 
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c 


c 

c — 


c 

c — 


c 

c — 


c 

c 

c 


c 

c 

c 

c 

c 

100 

c 


c 

3 

c 

c 


tal = (alf*t + alf *bta*t**2) / zO 

tml_p = 0. 

tml_t = (alf + 2 . *alf *bta*t ) / zO 

tm2 = p*pi/tau*phid(ttc) / Z 0 

tm2_p = pi/tau*phid(ttc) / Z 0 

tm2_t = p*pi/tau*phidd(ttc)*ttc_t / zO 

— residual 2: caloric equation 

res2 = h - (tml + tm2) 

r2_p = - (tml.p + tm2_p) 

r2_t = - (tml_t + tm2_t) 

— set Jacobian matrix 
a(l,l) = rl_t 

a(l ,2) = rl_p 

a(2,l) = r2_t 

a(2,2) = r2_p 

■” find inverse Jacobian matrix 

dotinv = 1.0 / (a(l ,l)*a(2,2) - a(l ,2)*a(2, 1)) 

ai(l ,1) = a(2,2)*detinv 

ai(2,2) = a(l,l)*datinv 

ai(l, 2) = -a(l ,2)*datinv 

ai(2,l) = -a(2 , l)*datinv 

~ set Hevton changes 
dt = -(aid ,l)*real + ai(l, 2 )*raa 2 ) 
dp = -(ai(2,l)*r«al + ai(2,2)*res2) 

rlx = 1.0 

if(rlx*dp .gt. 2.5*p) rlx = 2.5*p/dp 
if(rlx*dp .It. -.8*p) rlx = -.8*p/dp 
if(rlx*dt .gt. 2. Sat) rlx = 2.5*t/dt 
if(rlx*dt .It. -.8*t) rlx = -.8*t/dt 

' update variables 
t = t + rlx*dt 
p - p + rlx*dp 

■ convergence check 

if (aba (dp/p) .1*. apa .and. aba(dt/t) .la. apa) goto 3 
continue 

write(*,*) ’IGASPT: Convergence failed.’ 
vrita(«,«) 'dp dT dp, dt 
writa(*,*) 'p T h r: * , p, t, h, r 

continue 

set residual derivatives wrt input r,h variables 
rl.r = -p/ (r**2*t) 
rl_h = 0. 
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r2_r = 

0 . 

r2_h = 

1 . 

b(l,l) 

= rl_r 

b(l ,2) 

= rl_h 

b(2, 1) 

= r2_r 

b(2 ,2) 

= r2_h 


c set p,t derivatives wrt r,h 

t_r = -(ai(l,l)*b(l,l) + ai(l ,2)*b(2,l)) 
t_h = -(aid ,l)*b(l ,2) + ai(l ,2)*b(2,2)) 
p_r = ~(ai(2 , l)*b(l , 1) + ai(2,2)*b(2,l)) 
P-h = -(ai(2 , l)*b(l ,2) + ai(2,2)*b(2,2)) 


c 

c 


c 


c 


c 


c 


set second residual derivatives wrt r,h 
ttc = l./(tau*t) 
ttc_t = -1 ./(tau*t**2) 
ttc.tt = 2 . / (tau*t**3) 


z = 1. + p*pi*phi(ttc) 

z_p = pi*phi(ttc) 

z.pt = pi*phid(ttc)*ttc_t 

z_pp = 0. 

z.t = p*pi*phid(ttc)*ttc_t 

z_tt = p*pi*(phidd(ttc)*ttc_t**2 + phid(ttc)*ttc_tt) 


rl = p/(r*t) - z /z0 

rl_p = l./(r*t) - z _p /zO 

rl_pt = -l./(r*t**2) - z_pt/zO 

rl -PP = - z_pp/zO 

rl_t = -p/(r*t**2) - z_t /zO 

rl_tt = 2.*p/(r*t**3) - z_tt/zO 
rl_r » -p/(r**2*t) 

rl_h = 0. 
rl_hp = 0. 

rl_ht = o. 

rl.rp = -l./(r**2*t) 

rl_rt = p/(r**2*t**2) 

rl_rr = 2.*p/(r«*3*t) 


t«l = (alf*t + all *bta*t»*2) / zO 

t»l_t = (all + 2 . *alf *bta«t ) / zO 

* ( 2.*alf*bta ) / zO 

tal_pt S 0. 

tal_p = o. 

tal.pp a o. 


tn2 = p*pi/tau*pliid(ttc) / Z 0 
tm2_p * pi/tau*phid(ttc) / zO 
tm 2 _pt = pi/tan*phidd(ttc)*ttc_t / *0 
ta2_pp a 0. 

ta2_t = papi/tan* phidd(ttc)*ttc_t / zO 


tm 2 _tt * p*pi/tau*(pbiddd(ttc)*ttc_t **2 ♦ 
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* phidd(ttc)*ttc_tt) / zO 

c 

r2 = h - (tal + ta2) 

r2_p = - (tal.p + tm2_p) 

r2_t = - (tml_t + tm2_t ) 

r2_h = 1. 
c 
c 


c 


c 


c 

c 


c 


c 


c 



c 


- set and linearize new residuals: rlh = drl/dh = 0, r2h = dr2/dh 
ph = p_h 
th = t_h 

rlh = rl_p *ph + rl_t *th + rl.h 

rlh_ph = rl_p 

rlh_th = r 1 _t 

rlh.p = rl_pp*ph + rl_pt*th + rl_hp 

rlh.t = rl_pt*ph + rl_tt*th + rl_ht 

rlh_h = 0. 

rlh.r = -ph/(r**2*t) + th+p/(r**2*t**2) 


1. - tal_t*th - tal_p*ph - ta2_t*th - ta2_p*ph 
“ tml_p - tm2.p 

- tml.t - ta2_t 

- tal_pt*th - tml_ppeph - ta2_pt*th - tm2_pp*ph 

- tal_tt*th - tnl_pt*ph - tm2_tt*th - ta2 pt*ph 

0. 

0. 

rlh_th 
rlh.ph 
r2h_th 
r2h_ph 

8 1.0 / (a(l , l)*a(2,2) - a(l ,2)ea(2,l)) 

aih(l,l) = a(2 ,2)*detinv 

aih(2,2) = a(l , l)*detinv 
aih(l,2) = -a(l ,2)*detinv 
aih(2,l) - -a(2,l)*detinv 

dth = -(aih(l,l)*rlh + aih(l,2)*r2h) 
dph * -(aih(2 f l)*rlh ♦ aih(2,2)*r2h) 

ph ■ ph + dph 
th * th + dth 


r2h 

r2h_ph = 
r2h_th = 
r2h_p = 
r2h_t = 
r2h_h = 
r2h_r = 

a(l,l) = 
•( 1 , 2 ) = 
a(2,l) = 

a(2,2) = 

detinv 


set end linearize new residuals: rlr = drl/dr * 0, r2r = dr2/dr = 
pr = p.r 
tr = t.r 


rlr = rl.p *pr + rl.t *tr + rl.r 

rlr_pr * rl.p 

rlr.tr = rl.t 


= 0 


0 
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c 

c 


c 


c 


c 


c 

c 

c 

c 

c 


c 


c 


c 

c 


c 



c 


rlr_p = rl_pp*pr + rl_pt*tr + rl_rp 
Tix^t = rl_pt*pr + rl_tt*tr + rl_rt 
rlr - r = rl_rp*pr + rl_rt*tr + rl_rr 
rlr_h = 0. 


t2t = - tml_t *tr 
r2r_pr = 

r2r_tr = - tml.t 
r2r_p = - tml_pt*tr 
r2r_t = - tal_tt*tr 
r2r_h = 0. 
r2r_r = 0. 


tml_p *pr - ta2_t *tr - ta2_p *pr 
tml_p - ta2_p 

- tm2_t 

tml_pp*pr - tm2_pt*tr - tm2_pp*pr 
tal_pt*pr - ta2_tt*tr - tm2_pt*pr 


*(1,1) = rlr_tr 

a(l ,2) = rlr_pr 

*(2,1) = r2r_tr 

*(2,2) = r2r_pr 

dot in v = 1.0 / (a(l ,l)*a(2,2) - a(l,2)*a(2,l)) 

air(l,l) = a(2,2)*dotinv 

air(2 ,2) = a(l ,l)*detinv 

air(l ,2) = -a(l ,2)*dotinv 

air(2,l) = -a(2 , l)*detinv 


dtr = -(air(l,l)*rlr + air(l ,2)*r2r) 
dpr = -(air(2,l)*rlr + air(2,2)*r2r) 


pr = pr + dpr 
tr = tr + dtr 


calculate roaponaoa in dt/dh and dp/dh to unit h perturbation 
drlh = rlh_h + rlh_p*ph + rlh_t*th 
dr2h = r2h_h + r2h_p*ph + r2h_t*th 

drlr = rir_h + rlr_p*ph + rlr_t*th 
dr2r = r2r_h + r2r_p*ph + r2r_t*th 


dth = -(aih(l,l)*drlh + aih(l ,2)*dr2h) 
dpt * -(aih(2 , l)*drlh + aih(2,2)*dr2h) 
thb • dth 
phh * dph 

dth a -(air(l,l)*drlr + air(l,2)*dr2r) 
dph * ~(air(2 t l)*drlr + air(2,2)*dr2r) 
thr = dth 
phr = dph 


calculate responses in dt/dh and dp/dh to unit r perturbation 
drlh = rlh_r + rlh_p*pr + rlh_t*tr 
dx2h = r2h_r + r2h_p*pr + r2h_t*tr 
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drlr = rlr_r + rlr_p*pr + rlr_t*tr 
dr2r = r2r_r + r2r_p*pr + r2r_t*tr 
c 
c 

dth = -(aih(l ,l)*drlh + aih(l ,2)*dr2h) 
dph = -(aih(2,l)*drlh + aih(2 ,2)*dr2h) 
trh = dth 
prh = dph 
c 

dth = -(air(l,l)*drlr + air(l ,2)*dr2r) 
dph = -(air(2,l)*drlr + air(2,2)*dr2r) 
trr = dth 
prr = dph 

■at final first and second derivatives art (r,h) 

P- r = pr 

t_r = tr 

p_h = ph 

t_h = th 

p_hh = phh 

t_hh = thh 

p_rr = prr 

t_rr = trr 

p_rh = . S*(prh+phr) 

t_rh * .S*(trh+thx) 

return 

end 


subroutine nonstag(hO,rho,q, pO,pO_r.pO_q, 
* rO,rO_r,rO_q ) 


c Calculates stagnation pressure and density for 
c specified stagnation enthalpy, density, end speed. 


c 

Input : 


c 

hO 

stagnation anthalpy 

c 

rho 

danaity 

c 

c 

q 

spa ad 

c 

Output: 


c 

po 

stagnation pros sura 

c 

P0.r 

dpO/dr 

c 

po.q 

dpO/dq 

c 

rO 

stagnation dansity 

c 

rO_r 

dxO/dr 

c 

rO.q 

drO/dq 


c 


implicit real*4 (a-h,n,o-x) 
dinension a(2,2), ai(2,2), b(2,2) 
real*4 h_p,h_t 


3 ? 
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c 

c 


c 

ccc 

c 

c 

c 

c 

c 


c 


c 


c 


c 

c 

cc 

cc 

cc 

cc 

cc 

cc 

c 


comuon /noogai/ 

ft alf, bta, pi, tau, zO 

common /non* it/ 
ft c2, cl, cO 

data eps /5.0E-6/ 

z(pp,tt) = 1. + pp*pi*phi (1 ./(tau*tt)) 

Z-P^PP.tt) = pi*phi (1 ./(tau*tt)) 

z_t(pp,tt) = pp*pi*phid(l ./(tau*tt)) / (-tau*tt**2) 

h = hO - .5*q**2 

h_q = - q 

h_hO =1.0 

r = rho 

- sat input pressure and temperature and derivatives 
call ngaspt (h , r , p , p_r , p_h , p_rr , p_hh , p_rh , 

* t ,t_r,t_h,t_rr,t_hh,t_rh) 

- set entropy s and derivatives wrt p,t 

ttc = l./(tau*t) 

ttc_t = -1 ./(tau*t**2) 
ttc_tt = 2./(tau*t**3) 

ph = phi(ttc) 
phd = phid(ttc) 
phdd = phidd(ttc) 
phddd = phiddd(ttc) 

ph_t = phd * ttc_t 

phd_t = phdd * ttc_t 

phdd.t = phddd * ttc.t 

■ = alf*log(t) + 2.0*alf*bta*t 

k - p*pi*( t*phd *ttc_t + ph ) - log(p) 

s_p = - pi*( t*phd *ttc_t + ph ) - 1.0/p 

s_t = alf/t + 2.0*alf *bta 

k - p*pi*( phd *ttc_t + ph_t 
k + t *phd_t*ttc_t 

k + t*phd *ttc_tt ) 

' initial guess for pO,tO froa inperfect gas 
if (bta. eq. 0.0) then 
to = hO/alf 

else 

tO = (-1.0 + sqrt (1.0 + 4.0*bta*h0/alf)) / (2.0*bta) 
endif 

pO = p * erp(-alf*log(t) + alf *2.0*bta*(l .0-t) ) 
to = t 

pO = p 
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— lavton loop to converge on correct pO.tO 
itcon = 15 

do 100 iter=l, itcon 


ttc = l./(tau*t0) 

ttc_tO = -1 . /(tau*t0**2) 
ttc_ttO = 2 . /(tau*t0**3) 

ph = phi(ttc) 

phd = phid(ttc) 

phdd = phidd(ttc) 

phddd = phiddd(ttc) 

ph.tO = phd * ttc.tO 

phd.tO = phdd * ttc_tO 

phdd_tO = phddd * ttc_t0 

- enthalpy residual 

real - (alf*(tO + bta*t0**2) + pO*pi/tau*phd )/z0 

rl - p0 = ^ pi/tau*phd )/z0 

rl_t 0 = (alf*(l. 0 + bta*t0*2.) + pO*pi/tau*phd_tO)/zO 

- entropy residual 

rea2 = ali*log(tO) + 2 . 0*alf *bta*tO 
k - pO*pi*( tO*phd *ttc_tO + ph ) - log(pO) 

r2_p0 = - pi*( tOephd *ttc_tO + ph ) - 1.0/p0 

r2_t0 = alf/tO + 2.0*alf*bta 

* - p0*pi*( phd *ttc_tO + ph_tO 

k + tO*phd_tO*ttc_tO 

k + tO*phd *ttc_ttO ) 


- setup and invert Jacobian matrix 
&(1 » 1) = rl.tO 
*(1,2) = rl_pO 
*(2,1) = r2_t0 
*(2,2) = r2_p0 

detinv = 1.0 / (a(l,l)*a<2,2) - «(1 „2)*a(2, 1) ) 

aiCl.l) = »(2,2)*detinv 

ai(2,2) = a(l ,l)*detinv 

al(l, 2) = -a(l,2)*detinv 

*1(2,1) » -a(2,l)*detinv 

sat larton variables 

dt * -(al(l,l)*raal ♦ ai(l,2)*rea2) 
dp = -(al(2,l)*raal + ai(2,2)*re«2) 

rlx « l.o 

if(rlx*dp .gt. 2.5*p0) rlx * 2.S*pO/dp 
if (rlx*dp .It. - . 8*pO) rlx = -.8*p0/dp 
if(rlx*dt .gt. 2.5*t0) rlx » 2.5*t0/dt 
if(rlx*dt .It. ~.8*t0) rlx a -,8*t0/dt 


update variables 



pO = pO + rlx*dp 
tO = to + rlx*dt 
c 

c convergence check 

if (abs(dp/pO) .le. eps .and. abs(dt/tO) .le. epa) go to 2 
c 

100 continue 

c 

write(*,*) ’H0SST1G: Convergence failure.’ 
write(* , *) ’dp dT :’,dp, dt 
write(* ,*) >po To h r: ’ ,pO,tO,h,r 
c 

2 continue 

c 

c set residual derivatives urt (s,hO) 

rl_a = 0. 

r2_s = -1.0 

rl_h = -1.0 

r2_h =0. 
c 

b(l,l) = rl_a 

b(l ,2) = rl_h 

b(2,l) = r2_a 

b(2,2) = r2_h 

c 

c set (t0,p0) derivatives wrt (s,hO) 

tO.s = -(ai(l ,l)*b(l , 1 ) + ai(l,2)*b(2,l)) 

ccc to.ho = -(ai(l , l)*b(l , 2 ) + ai(l , 2 )*b( 2 , 2 ) ) 

pO_s = -(ai(2,l)*b(l,l) + ai(2,2)*b(2,l)) 

ccc pO.hO = -(ai(2,l)*b(l,2) + ai(2,2)*b(2,2)) 

c 

c convert derivatives wrt (s.hO) to wrt (p.t.hO) 

to_t = to_**§_t 

tO_p = tO_**g_p 
pO_t = pO_i*»_t 
pO_p = pO_**e_p 

c 

c 

c set stagnation density rO and derivatives wrt (pO.tO) 
zz = z(pO,tO) 

«-P = z_p(pO,tO) 
zz_t * z_t(pO,tO) 
c 

rO * zO/zz * pO/tO 

rO_z = -z0/zz**2 * pO/tO 

c 

rO_pO = rO_z*zz_p + zO/(zz*tO) 

rO_tO = rO_z*zz_t - zO*pO/(zz*tO**2) 

c 

c convert derivatives fron wrt (pO.tO) to wrt (p.t.hO) 

rO_p a rO_pO*pO_p + rO_tO*tO_p 

rO_t = rO_pO*pO_t + rO_tO*tO_t 

ccc rO_hO = rO_pO*pO_hO + rO_tO*tO_hO 
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c convert derivative* from wrt (p,t) to wrt (r,q,hO) 
rO.r = rO_p*p_ r + rO_t*t_r 
rO_q = (rO_p*p_h + rO_t*t_h) *h_q 
pO_r = pO_p*p_r + pO_t *t_r 
pO_q = (pO_p*p_h + pO_t*t_h)*h_q 
c 

ccc rO.hO = (rO_p*p_h + rO_t*t_h) *h_hO + rO^hO 

ccc pO_hO = (pO_p*p_h + pO_t*t_h) *h_hO + pO.hO 

c 

return 

end 


c- 

c 

c 

c- 


real*4 function phi(ttc) 
implicit real*4(a-h,m,o-z) 


Return* function phi used in non~ideality parameter 
Z = 1 + pi*phi(ttc) 

common /nonfit/ 

* c2 , cl, cO 


phi = c2*ttc**2 + cl*ttc + cO 

return 

end 


real*4 function phid(ttc) 
implicit real*4(a-h t m,o-z) 
common /nonfit/ 

* c2, cl, cO 

phid = 2 . *c2*ttc + cl 

return 

end 


real*4 function phidd(ttc) 
implicit real*4(a-h,m,o-z) 
common /nonfit/ 

^ c2 f cl , cO 

phidd * 2 . *c2 

return 

end 


real*4 function phiddd(ttc) 
implicit real*4(a-h,m,o-z) 
common /nonfit/ 

^ c2 , cl i cO 
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phiddd = 0. 

return 

end 


subroutine hgent (hO ,r ,q, s) 


c Returns entropy s from input variables hO,r,q 



common /nongas/ 

* alf , bta, pi, tau, zO 
common /nonlit/ 

* c2, cl, cO 
c 

h = hO - . 5*q**2 
c 

c set input pressure and temperature and derivatives 
call ngaspt (h , r , p , p_r , p_h , p_rr , p_hh , p_rh , 

* * , t_r , t_h , t _rr , t_hh , t_rh) 
c 

ttc = l./(tau*t) 
ttc_t = -1 . /(tau*t**2) 
c 

ph = phi (ttc) 

phd = phid(ttc) 

c 

s = alf *log(t ) + 2.0*alf*bta*t 

* - p*pi*(t*phd*ttc_t + ph) - log(p) 

return 

end 


subroutine nonga*v(h0 ,r,q, gaa.gan.r ,gaa_q) 



c Returns •'equivalent" game for BL density profile 

_____ 

common /nongas/ 

* all, bta, pi, tau, zO 
common /nonfit/ 

* c2, cl, cO 
c 

c set B tatic enthalpy 

h = hO - O.S*q*e2 
h.q = -q 

c 

c set pressure and temperature and derivatives 

call ngaspt(h.r,p,p_r,p_h,p_rr,p_hh,p_rh, 

* t .t_r,t_h,t_rr ,t_hh,t_rh) 
c 

c set speed of sound squared: a*2 = dp/dr (at constant s) 

»«q ■ p_r / (1. - p_h/r) 
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•■q^r = p_rr / (1. - p_h/r) 

* ~ P_r / (1. - p_h/r)**2 *(p_h/r**2 - p rh/r) 

•■q_h = p_rh / (1. - p_h/r) 

* + P_r / (l. - p_h/r)**2 *p_hh/r 

c 

c 

ttc = l./(tau*t) 
ttc.t = -1 . /(tau*t**2) 

= 2 . /(tau*t**3) 

c 

ph = phi(ttc) 
phd = phid(ttc) 
phdd = phidd(ttc) 
phddd - phiddd(ttc) 


c 


c 


c 

c 


c 


c 


z = 1 . + p*pi*ph 
Z -P ~ pi*ph 

z -t = p*pi*phd*ttc_t 

cp = ( alf *( 1 . 0 + 2.0*bta*t) 

* + p*pi/tau* phdd*ttc_t ) / Z 0 

C P-P = ( pi/tau* phdd*ttc_t ) / zO 
cp_t = ( alf *( 2.0*bta ) 

* + p*pi/tau*(phddd*ttc_t**2 + phdd*ttc_tt) ) / zO 

Z6t = h/(cp*t)*(l .0 - p*pi/ (t*tau)*phd/z) * zO 

zat.h = 1.0/(cp*t)*(l.o - p*pi/(t*tau)*plid/z) • zO 

Z#t -P = h /(c P *t)*( - pi/ ( t *t au) *phd/z 

* , , % - P*pi/(t*tau)*plid/z*<-z_p/z)) * z0 

* ~ (zat/cp) *cp_p 

zat_t = h/(cp*t)*( - p* P i/(t*t. a )*phd/z*(-z_t/z - i.o/t) 

- p*pi/ (t*tau)*phdd*ttc t/z ) * zO 

* - (zat/cp) *cp_t - (zat/t ) 


g“ = asq/(h*zat) +1.0 
g“-r = 

gaa_h = aaq/(h*zat)*(-zat_h/zat - 1.0/h) 
g««_p = asq/ (h*zat ) * (-zat.p/zat ) 
g“-t = *sq/(h*zat)*(-zat_t/zat) 


■«q_r/(h*zat) 
+ asq_h/(h*zat) 


g«M_h * g**»_p*p_h + gu.t»t.h + gaa_h 
*“- r 3 ! tt -P*P- r ♦ gan_t*t_r + gaa.r 


g“-q * gan_h*h_q 

return 

end 


snbroutina sonic (hO.pO.rO, q,p,r) 

calculator sonic qnantltlas q,p,r 
fro* spacifiad sonic qoantitlas hO.pO.rO 
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implicit rami (m) 
data aps / l.Oa-6 / 

c 

c initialize with perfect gas 

gam = r0*h0 / (r0*h0 - pO) 
gml = gam -1.0 
c 

q = sqrt (2 . 0*h0/ (2 . 0/gml + 1.0)) 


c 

c — 


c 


c 


c 


c 


trat = 1.0 + 0 . 5*gml 
p = pO*trat**(-gam/gml) 
r = r0*trat**(-l .0/gml) 


converge on non-ideal values by forcing H~2 = 1, and pstag = pO 
do 10 iters=l, 15 

call nideal(h0 ,r,q, p ,p_r ,p_q, 

* msq,msq_r ,msq_q ) 
call nonstag(hO f r,q, pstag ( pstag_r,patag_q, 

* rstag,rstag_r ,rstag_q ) 
real = msq - 1.0 

all = msq_r 
al2 = msq_q 

re§2 = pstag - pO 
a21 = pstag.r 

a22 = pstag_q 

detinv = 1.0/(all*a22 - al2*a21) 
dr = -(resl*a22 - al2 *re*2)*detinv 
dq = -(all *res2 - resl*a21 )*detinv 

dp = p_r*dx + p_q*dq 

rlx = 1.0 

if (rlx*dr .gt. 1.5*r) rlx = 1.5*r/dx 
if (rlx*dr .It. -.e*r) rlx = -.6*r/dr 
if (rlx*dq .gt. 1.5*q) rlx =* l.$*q/dq 
if (rlx*dq .It. -,6*q) rlx = -.6*q/dq 


r = r + rlx*dr 
q 3 q + rlx*dq 
p * p + rlx*dp 

<ta« = amaxl ( abs(dr)/r , aba(dq)/q ) 
if (dmax .It. apa) go to 11 

10 continue 

write(* f *) * sonic : convergence failed, dmax =>, d max 

11 continue 


return 
end ! sonic 
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Appendix B 

High-Order Airfoil Farfield Boundary Conditions 
for Ideal and Non-Ideal Gas Flows 


The steady flow around an airfoil away from shock wakes and viscous regions has constant 
entropy and total enthalpy, and hence is also irrotational. These properties hold whether the fluid 
is an ideal or a non-ideal gas. The flow can then still be decribed by the velocity potential $ or 
the perturbation potential <fi. Assuming the freestream is aligned with the z-axis, the following 
relations are obtained. 


V<£ = q ~ 

<? 2 = k1 2 = 

iv(g 2 ) = qVq = 


.(* + <t>) 

Qoo [( 1 “I" 4>x)i F &yj] 

qL [(i + <t>x) 2 + 4>l\ 

qll (4> xx F <t>x <f> XX F <t>y<t>xy)i 
F {$xy F tyx <f>xy F <t>y<f>yy)j } 


The governing flow equation is: 


V-(pV$) 
or V 2 $ 



In isentropic flow (5 


and hence 


constant), p = p(p), so 
dp 


Vp = 


dp 


Vp 


-4 1*1 


V 2 $ 
a 2 V 2 (j> 


4 qVq ■ 
a L 

qVq • [(1 + 4 > x )i + <f> y j ] 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 


( 7 ) 


( 8 ) 

( 9 ) 


where a is the speed of sound. In isentropic, adiabatic flow, the speed of sound is uniquely related 
to the speed: a = a(q). For a perfect gas, a(q) is given by 

do) 

while for an imperfect and/or non-ideal gas it is necessary here to linearize a(q) about the freestream 
conditions. 


2 2 

a 2 ~ ai + v °° ; 


d(ll) 




It is convenient to define an “equivalent” ratio of specific heats 7' for the non-ideal gas as 

,<*(<£)! 


7 ' = 1 


d(q i) 


( 11 ) 


( 12 ) 
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so that the a{q) relation for the non-ideal gas can be compactly written as 


«’ = <>l - ^ - <d) (13) 

For a perfect gas, this reverts to the exact form (10) since in this case 7' = 7. It is interesting 
to note that 7' can easily be less than unity for heavy gases such as sulfur hexafluoride, while 
invariably 7 > 1 for perfect gases. 

Substituting for a 2 , <7 2 , and qVq in equation (9), we obtain 



o * 00 


( 20* + 4>l + 


[<t>xx + <f>: 


vvi 


Qoo [(^ 4 4*x) {$xx 4 <^x ^11 4" 0y<^xy) 


4 ^y (0xy 4 <^x <^xy 4 ^y^yy)] (14) 



[^X 4 <£yy] 


[<£ xx + 2<f> x <t> xx + 2<j>y<j> xy \ + O 


(15) 


(1 4 [(7 I +1) < / > X0XX 4 ( 7 1 )<Ax <£yy 4 24 > y <f> xy ] + o(f) (16) 

where M = q^fa^ is the freestream Mach number. 

Equation (16) is the 2D second- order Prandtl-Glauert equation which governs small-perturbation 
non-ideal compressible potential flows. It has the same form as the equation for a perfect ideal gas 
as derived in references [4] and [5], except that the usual ratio of specific heats 7 is replaced by 
the “equivalent” value 7' defined by equation (12). Wagner and Schmidt [1] have considered the 
first-order version of equation (16) using 7' in lieu of 7. 

In terms of the Prandtl-Glauert coordinates 


x = x//3 

y 


y 

-2 


= * 2 + y 2 


9 — arc tan — 
x 


where (3 = >/l — M£, the general solution to equation (16) is 

<f> = 


-r e 

—0 + -^lnf 
2tc 2t 


4 

4 


D x cos 9 £L sin 9 


where 


k x 


2x r 

/ PMx 

V 2tt 

3-7' 


4 


2w T 


)*[■ 


2 r , lnr cos 30 
k\ — — «2 


I ( 7 ' +1 4. 3 ~7' \ 

4 V /3 3 5 ) 


/3 3 p 

Terms of order 1/f 2 and above have been discarded. 


J_ /r+i , 7'~1 \ 
16 V /3 3 /9 /' 


(17) 

(18) 

(19) 

(20) 


( 21 ) 

( 22 ) 
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In a flow solver, the circulation T can be determined either directly from the lift per unit span 
V (Euler or Navier-Stokes code), 

r = A 

PooQ oo 

or indirectly by specifying a Kutta condition (potential solver or MSES). The source strength E 
can be determined from the total profile drag per unit span D', or from the asymptotic mass defect 
behind the airfoil including the shock wake. 



D' 

P~vl 


(24) 


In the case of a potential solver, D f should not include the wave drag since there is no shock wake 
(unless an entropy correction scheme is employed). Note that T and E here have units of length 
since <j> in (1) corresponds to a unit freestream speed. 

Cole and Cook [5] give explicit expressions for the doublet coefficients D x and D y in terms of 
field integrals over the domain. Unfortunately, these expressions are unwieldy and for a non-ideal 
gas would be rather expensive. A simpler and economical approach is to iteratively update D x and 
D y by minimizing the mismatch between V# and the velocity Solution from the flow solver on the 
outer boundary. The approach taken in reference [4], for example, is to minimize the integral 


/ = 



X ^solution | da? 


(25) 


taken over the outermost streamlines. The doublet terms in the farfield expansion (21) decay 
faster than the others, and so can be neglected for sufficiently distant outer boundaries. However, 
retaining them greatly reduces the sensitivity of the solution to domain size, especially for transonic 
flows [6]. 

With its term coefficients defined, equation (21) gives a very accurate representation of the 
perturbation potential <t> away from the airfoil. The gradient of equation (21) accurately gives 
the total velocity q via relation (2). Either <^>, q \ or an appropriate derived quantity may then be 
imposed at the outer domain boundary as a high-order boundary condition. A potential solver 
would typically impose <j> or d<t>jdn , whereas an Euler or Navier-Stokes solver would typically 
impose the flow angle at the inflow find pressure at the outflow, both being determined from V$. 
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Appendix C 

Shape Parameter Relations 
for Ideal and Non-Ideal Gas Flows 


The major influence of compressibility on boundary layer behavior is a non-uniform density 
profile, which alters the layer’s response to pressure gradients. In an integral scheme for adiabatic 
flows, this effect is mostly captured by the correlation between shape parameter H , the kinematic 
shape parameter H k, and the edge Mach number M e . The shape parameters are defined as 

H = I(l~RU)dy = f(l-U)dy 

J(l-U)RUdy k J(l~U)Udy 

where U(y) and R(y) are the velocity and density profiles. 




Ue 



(27) 


Since the velocity profile U(y) and hence Hk are only weakly affected by compressibility, reduction 
of the density profile R(y) near the wall due to adiabatic heating will increase H as can be seen 
from its definition (26). In turn, the von-Karman integral momentum equation 


dQ _ C f 
dx 2 


(2 + H -M e 2 ) 


0 du e 
u e dx 


(28) 


shows that an increase in H will increase the momentum thickness growth rate dO/dx for a given 
adverse pressure gradient. The integral boundary layer formulation in MSES (and its precursor 
ISES [4]) employs a correlation of the form Hk{H, M e ) for air. This is re-derived for the non-ideal 
gas model as follows. 


As developed in Appendix A, the state equation of a non-ideal gas can be written as 


P 

pTZT 


Z(Pi t) 


while the corresponding caloric equation in differential form is 


(29) 


dk = c p {r)dT + d\pT{r )] - c p {r)dT + pT\r)dT = c p (p,r)dT (30) 

where the approximation is made on the basis that dp ~ 0 across a boundary layer. Linearizing 
the caloric equation across the boundary layer we have 


h-K = c p (T-T e ) 

T \h ) c p T 

(31) 

(K \ he 


\h V c Pt T e 

(32) 

The non-idecility factor Z for most non-ideal gases has the form 


Z( p ,t) = l + —4>(t c /t) 

Pc 

(33) 


3t 


n 



where p c and T c are the critical pressure and temperature. This can likewise be linearized about 
the edge conditions as follows. 


Z - Z' 

pe jt f T c T c \ 

(34) 

Z' 

z 

Pc T c (T t \ 4>’ e 

Pc T'\T ) Z 

(35) 

Combining this with equation (32), we 

have 


*-> 

P« e T c fh t \ h e # 

(36) 

Pc T e \ h ) c Pt T e Z e 


Using the equation of state (29), the density profile is then related to the T and Z profiles as 


P_ _ 

Pe ~ T Z Pe 
T Z 


(37) 


with the usual boundary layer approximation p ^ p c being made. Using relations (32) and (36), 
the density profile can be written in terms of the enthalpy profile alone. 


where 


R = 

1 + 

(*- 

0 

1 Pe T c (K A K <t>' e 

Pc T e \ h ) c Pt T e Z e 

(38) 

= 1 + ( 


) « ‘n (> 

-7'W 

2 ] (39) 

R 1+ ( 

T ~ 1 

)‘ 


(40) 



C 

‘ Cpje 1 

Pe ^ K \ 
\ PcT'Zj 

(41) 


For turbulent adiabatic boundary layer flows, it is reasonable to assume a constant stagnation 
enthalpy across the layer, although this is strictly true only for a turbulent Prandtl number of 
unity. Since the turbulent diffusion mechanisms of momentum and heat are essentially the same 
in a gas (convection by eddies), turbulent Prandtl numbers are typically close to unity. Hence, 
the assumption of constant stagnation enthalpy is reasonable. With h 0 denoting the stagnation 
enthalpy, the velocity and static enthalpy profiles are then related by 


he _ ho - u 2 J2 _ 1 - 

h h 0 - u 2 /2 i _ JjL (J 2 

**_! = 
h 

and the density and velocity profiles are then related by 


R = 1 + 






(u 2 - i)C- 


(42) 

(43) 


(44) 
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Since u\jh o and ( are both functions of the edge Mach number M c , the density profile (44) 
implicitly defines H * in terms of H and M e . To obtain this relation in closed form, it is necessary 
to assume a small- defect profile 


U = 1 - c ; e < 1 

so that the density profile can be approximated by 

R = 1 + — %-(-2<r)C + 0(e 2 ) 

1 ^ 

R = 1 - (l/v-l)M 2 e( 

with the convenient “viscous” equivalent ratio of specific heats *) v defined by 

^M 2 ( = 


1 - 




= 1 + 


a 2 1 


h 0 - u\j 2 ( 

with a e being the speed of sound at the boundary layer edge. 

The shape parameter H now becomes 

H = /{l-[l-(7»-l)MM(l-6)} dy 

J{e[l-( lv -l)M 2 e](l-e)}dy 


Jedy + (7„-l)M e 2 /€(l 

Ml-‘)dy - ( lv -l)M?Jt2(l-e)dy 


Jedy 


+ ( T „-l )M 2 + 0(e 2 ) 


fe(l - e) dy 
* H k + ( 7 „-l)M e 2 

The required shape parameter correlation is therefore 

H k = H - ( 7 „-l )M 2 . 


(45) 

(46) 

(47) 

(48) 

(49) 


(50) 

(51) 


In the limiting case of a perfect gas, 7„ = 7. For 7 = 1.4 (air), MSES presently uses Whitfield’s 
correlation [7] in this case is 


H - 0.29 M 2 
1 + 0.113M* 

= H - 0AM 2 + 0(M *) 


(52) 

(53) 


which is seen to be consistent with the more general non-ideal gas result (51). Whitfield’s particular 
form (52), however, is reportedly more accurate for Prandtl numbers somewhat less them unity 
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where the total enthalpy profile is not quite uniform as was assumed here. It is therefore appropriate 
to put correlation (51) into Whitfield’s form, while also incorporating the Prandtl number. The 
final shape parameter correlation is 


Hk{H,M e ) 


H - Pr(y v -l)M' 

1 + (1-Pr)( 7 „-1)M* 


(54) 


which reduces to Whitfield’s form for 7„ = 1.4, £ = 1, and Pr = 0.7, and to the non-ideal gas form 
(51) for Pr = 1 which was assumed in its derivation. 


It noteworthy that for most heavy gases 7„-l is considerably smaller them for air. For SFq with 
stagnation conditions at STP and M e = 1, for example, 7„-l = 0.17 for SF 6 and 7„-l = 7-I = 0.4 
for air. Hence, the influence of M e in SF 6 is smaller, and H values near a shock in SF 6 will be 
smaller them those in eiir. The smaller H values in turn reduce the boundary layer’s response to 
adverse pressure gradients as discussed above. The airfoil will therefore be more resistemt to Mach 
drag-divergence in SF 6 them in air. 


For simplicity, the implementation of the shape parameter correlation (54) in MSES assumes 
that 7„ is constemt, being evaluated from (49) at sonic edge conditions: a t = u e = a* , M t = 1. 
Given the degree of approximation used in deriving (54), it is felt that neglecting the edready weak 
dependence of ~f v on u e is justified. Freezing ~j v at the sonic conditions is judged appropriate since 
its effect on H becomes significant only for M e close to unity. 
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